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Among the sources of systematic error in nucleon structure calculations is contamination from 
unwanted excited states. In order to measure this systematic error, we vary the operator insertion 
time and source-sink separation independently. We compute observables for three source-sink 
separations between 0.93 fm and 1 .39 fm using clover-improved Wilson fermions and pion masses 
as low as 150 MeV. We explore the use of a two-state model fit to subtract off the contribution 
from excited states. 
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1. Introduction 

As Lattice QCD calculations of nucleon structure have been performed with decreasing pion 
mass, chiral extrapolations have become less able to reconcile lattice data with experiment [1,2]. 
An important source of systematic error is contamination from excited states [3, 4], i.e. the failure to 
isolate the ground-state nucleon. 

In order to compute nucleon matrix elements and form-factors, we compute nucleon two-point 
functions C2pt{T,p') and three-point functions Cf^^{T,z,p,p') on the lattice [5], where T is the 
Euclidean time separation between the source and the sink, T is the Euclidean time separation 
between the source and the operator ff, p' is the sink momentum, and p is the source momentum. If 
the nucleon interpolating operator with momentum p creates A'^ states with energies En{p), then the 
two-point and three-point functions have the form 

C2AT,p)=Z^n{p)^-^"^'^^ 

n,n'=0 i 

where t = —{p' — p)^, F"^"' {t) are transition form-factors for the transition from state n to state n' 
via the operator ^, and M is determined by kinematics and by the definition of the form-factors. 



2. Computing form-factors 

The traditional approach for computing nucleon matrix elements and form-factors is to use a 
ratio of three -point and two-point functions: 



^^^^ ^ ^ _ C%,{T,Z,p,p') jC2^,{T-T,p)C2,,{T,p') 



v/C2pt(r,p)C2pt(r,pO y C2p,(r - x,p')C2^,{x,p) ' 

such that in the limit of large T and T — z, where excited states have negligible contribution, the 
amplitudes aQ{p), ao{p') and the exponential dependences on z and 7 — T are canceled. In practice, 
we compute these ratios for one or more fixed values of T and T G [0, T], then for each (T, t) we 
compute the form-factors Fi{t) from an overdetermined fit to the ratios 7?. Then we produce "plateau 
plots" of Fi{t) versus T, for fixed T, and look for a flat region where it is hoped that excited-state 
contamination is small. 

The new approach introduced here is to do a combined fit to C2pt and C^pt using an A/^-state 
model derived from Eq. 1 . 1 , for a small value of A'^. In order to further simplify the fit, we assume 
the dispersion relation E^{p) = ml + p^, and define 

Fp^"'{t) = ^an{p)aAp')Fr"\t), 

so that for fixed masses, the fit model for C2pt and C^pt depends linearly on the fit parameters a„ {p) 
and F!^~^"' it). This enables the use of an iterative fitting procedure for the masses, where at each step 
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the larger set of linear fit parameters is solved for exactly. Denoting the two-point and three-point 
functions to which we are fitting by Ca and the linear fit parameters by bi, we minimize 



where 5*^ is our estimate of the covariance matrix of {Ca}- Because typically the number of 
independent samples that we have is of the same order as the number of variables Ca, we use 
a shrinkage estimator of the covariance matrix [6], S* = (1 —X)S + XT, where S is the sample 
covariance matrix, T is its diagonal part, and A is estimated from the data such that the expected 



If the number of states included in the model is less than the number of states that have 
non-negligible contributions to C2pt and C3pt, then each model state will have to account for the 
contributions from more than one state in the lattice data. The best-fit masses in the model will 
be weighted averages of the masses of states that contribute to the data. The relevant weights are 
different for different correlators [8], so when doing a 2-state fit, we use a single ground-state mass 

(2) 

niQ but allow the excited-state mass in the model to vary, using one mass m\ for C2pt and a different 

(3) 

mass m\ for C3pt. The difference between these two masses is an indicator of the importance of 
omitted states. 

3. Lattice measurements 

We use 2-1-1 flavors of tree-level clover-improved Wilson fermions coupled to double-HEX- 
smeared gauge links [9] . Results presented here are from ensembles with lattice spacing a = 0. 1 1 6 fm 
and pion masses approximately 150, 200, and 250 MeV. At the lightest pion mass, the lattice volume 
is 48^, and at the other two it is 32^ x 48. In order to study excited-state contamination, we compute 
nucleon three -point functions with three different source-sink separations T /a = S, 10, and 12. 
We focus on isovector quantities to avoid contributions from disconnected diagrams. We use the 
standard nucleon operator, Na = s"^'^ {u^Cy^dhlUca, where the quark fields have smearing tuned 
to minimize excited states seen in the two-point function, and we use the spin-parity projection 
Tpoi = 1-0^75 ^ p^w Qf j]^g errors below are statistical and are computed using the jackknife 
method. 

Matrix elements of the vector current Vq = qj^q are parametrized by the vector form-factors 
^1,2(0- The mean squared Dirac radius is defined from the slope of Fi{t) at zero t: F\{t) = 
F{{0)[\ — g(ri)^? + ^(?^)]. We compute the isovector Dirac radius from a linear fit to Fi(0) and 
F\{ti), where t\ corresponds to three -point functions with p equivalent to ^(1,0,0), and p' = 0. 
Using the ratio method, we average the three central points of each plateau plot to arrive at a value of 
{r"^'^Y foi" ^^'^h source-sink separation and each ensemble. The results are shown in Fig. 1. There 
is a consistent trend across the three ensembles of the Dirac radius increasing between T/a = 8 
and T /a = 10. However, the behavior going to T/a = 12 is ambiguous and it is unclear, from this 
analysis and with currently available statistics, whether or not excited-state contamination is still a 
problem at T/a = 10. 

We also compute the isovector average momentum fraction [x)"^'^ from forward matrix ele- 
ments of the operator Gq"^ = qy^^iD'^^q, where the braces denote taking the symmetric traceless 



= {Aai{{mr,])bi-CaW)-^l{Apj{{m,])bj-Cf,) 



error of the correlation matrix R^^ = S*^n/ ^ S*^a^*RR asymptotically minimized [7]. 
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Figure 1: Isovector Dkac radius (r"^ )^ (fm^) 
versus rrijc (GeV). At each pion mass, from left 
to right are results obtained from the fitting pro- 
cedure, and from the ratio method with three 
increasing source-sink separations. 
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Figure 2: Isovector average momentum fraction 
(x)"^'^ (bare) versus m^c (GeV). At each pion 
mass, from left to right are results obtained from 
the fitting procedure, and from the ratio method 
with three increasing source-sink separations. 



part of the tensor Note that the results presented here have not been renormalized. Again averaging 
the three central points of each plateau plot, the results are seen in Fig. 2. In this case, a clearer trend 
is visible: {x)"^^' decreases as T/a increases, and there is no evidence that even T/a = 12 might be 
large enough that excited-state contamination is negligible. This result is consistent with a recent 
study using the open sink method, where x is fixed in order to allow the calculation of three-point 
functions at all values of 7 [10]. Furthermore, the effect of excited states appears to be larger at 
smaller pion masses. This is consistent with the fact that earlier calculations at larger pion masses 
using domain wall valence fermions on an asqtad sea yielded chiral extrapolations to the physical 
point in excellent agreement with experiment [11, 12]. 



4. Fit results 



We perform 2-state fits to compute (r" ) . This involves fitting to C2pt and C^'^^'' , and the 
relevant (transition) form-factors are defined by [13]: 



{N„,{p',k')\V^\N„{p,k))=Unip',k') 



Un{p,X), 



where = p' — p. Combining 3-point functions that (up to an overall sign) have the same set of 
contributions from form-factors, at the first nonzero momentum transfer there are three independent 
^Spl (-^' which correspond to the three matrix elements listed in Fig. 3. Including these in a fit with 
the three-point functions at zero momentum transfer (both at rest and boosted), forT/a€ [\,T /a—\], 
T/a G {8, 10, 12}, and with C2pt(r,p), for p G {0, ^(1,0,0)} and T/a G [2, 12], we arrive at a fit 

(2) (3) 

to 211 variables, with 20 linear fit parameters and 3 mass parameters wiq, m\ , and m\ . As before, 
the Dirac radius is determined using a linear fit to the resulting F\ (0) and Fi [t]). 

Selected parameters and derived quantities from the fit to the niji = 250 MeV ensemble are 
shown in Tab. 1 . In particular, note that the best-fit excited-state mass for the three-point function is 
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zVdof 
A 

moa 

(2) 

m\ a 



a„ (0,0,0) 



,0,0) 



^r"'(fi) 



Fi{h) 
Fiih) 



89(22)/188 
0.036(12) 
0.637(6) 
1.55(14) 
0.99(5) 
2.41(9) 
^2.87(47)^ 

2.08(8) 
^2.49(41) 

2.09(9) 
-0.21(4) 

7.3(5) 
-0.5(3) - 
0.934(10) 
3.25(13) 



X 10 



-10 



X 10-1" 

-0.18(5) ' 
-0.23(27)^ 

-0.7(2) i 
11.1(7.4) i 



X 10 



-10 



X 10 



-10 



Table 1: Selected results from fit used to 

compute (r"-'')2 for the m„ = 250 MeV 

ensemble. Form-factors ^1,2(^1) are not 

renormalized and are computed from 
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Figure 3: Three-point function (points) and fit 
(error bands) versus z/a, with T /a = S. Matrix 
element labels are representatives from the sets 
of equivalent three-point functions that are aver- 
aged to compute the points shown here. Fit 

bands are determined by the quantities moa, 

(3) ~ ~ 
m\ a, F\{ti), and ^2(^2) listed in Tab. 1. Note 

that the points have correlated errors, and that 
neglecting correlations will cause the fit to over- 
lap with the data. 



lower than that for the two-point function. This suggests that there may be a state that has important 
contributions to three-point functions but is not easily detected from the two-point function of 
a single smeared nucleon operator. The fit model is compared with the 3-point functions at the 
first nonzero momentum transfer and T/a = 8 in Fig. 3. In this figure, contributions from the 
ground-state nucleon will decay approximately as exp(— 0.03T/a), since the nucleon with nonzero 
momentum at the source has slightly higher energy than the nucleon at rest at the sink. This is 
the approximate behavior seen in the first two of the three 3-point functions, however the third 
approaches zero more rapidly. This can be explained if the different three-point functions have 
different relative contributions from excited states. This is extra information that the fit is able to 
make use of, but is discarded when producing a plateau plot. 

In Fig. 1, we compare {r'l^'')^ from this fit with the values from the ratio method. The fit results 
show a stronger trend of increasing Dirac radius at smaller pion masses, although since the fit points 
have large errors (only slightly smaller than for the T/a = 12 plateau values), the outcome from the 
fit is consistent with with the T /a = 10 plateau values across all ensembles. 

To determine {x)"^'^, we perform the analogous fit to C2pt and C^^['', except that we restrict to 
r/a [2,T /a — 2], since the operator extends in the time direction. The (transition) matrix elements 
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are parametrized by five generalized form-factors: 



{NAp'^')W^'''^Nn{p,X))=uAp'^') 



+ 



2 



'20 



t) + Ai^^ : ^21 



Un{p,X), 



where p = {p' + p) /2, and A21 and B21 vanish for non-transition matrix elements due to time-reversal 
symmetry [14]. Including the first nonzero momentum tiansfer helps to give a better handle on 
the excited-state mass. Overall, the fit has 358 variables, with 38 linear fit parameters and the 
three masses. The fit results for (a:)"^'^ = A2o(0) are compared with the ratio method in Fig. 2. 
These show the same downward trend as the pion mass decreases that was seen in the r/<3 = 10,12 
plateau values, but with larger errors. At mji = 150 MeV, statistics were too poor, such that the value 
obtained from the fit, {x)"-'^ = -0.03(67), is not useful. 



5. Conclusion 



There is strong evidence of excited-state contamination in nucleon structure calculations. We 
see that the effect depends on the observable and consequently there is a clear need to quantify the 
systematic error from excited states for each observable. In particular, as demonstrated in Fig. 2, 
the use of a single source-sink separation smaller than about 1.4 fm is inadequate for accurate 
calculations of the isovector average momentum fraction near the physical pion mass. 

We have presented one promising approach to reducing this problem by including an additional 
finite number of excited states in a fit. Applying this method with one excited state, we obtain results 
that, with currently available statistics, are consistent with the ratio method. 
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